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Abstract. The generally held view that a model of large- 
scale structure, formed by collisionless matter in the Uni- 
verse, can be based on the matter model "dust" fails in the 
presence of multi-stream flow, i.e., velocity dispersion. We 
argue that models for large-scale structure should rather 
be constructed for a flow which describes the average mo- 
tion of a multi-stream system. We present a clearcut rea- 
soning how to approach the problem and derive an evolu- 
tion equation for the mean peculiar-velocity relative to 
background solutions of Friedmann-Lemaitre type. We 
consider restrictions of the nonlinear problem and show 
that the effect of velocity dispersion gives rise to an effec- 
tive viscosity of non-dissipative gravitational origin. We 
discuss subcases which arise naturally from this approach: 
the "sticky particle model" and the "adhesion approxi- 
mation" . We also construct a novel approximation that 
features adhesive action in the multi stream regime while 
conserving momentum, which was considered a drawback 
of the standard approximation based on Burgers' equa- 
tion. We finally argue that the assumptions made to ob- 
tain these models should be relaxed and we discuss how 
this can be achieved. 

Key words: Gravitation; Hydrodynamics; Instabilities; 
Methods: analytical; Cosmology: theory; large-scale struc- 
ture of Universe 



1. Why model a multi stream flow ? 

The problem of how to treat a multi-stream flow and 
the question of how to follow approximations beyond 
shell-crossing time are themes that are often repeated in 
the history of constructing models for large-scale struc- 
ture. The application of approximation schemes such as, 
e.g., Lagrangian perturbation solutions (which, to first or- 
der, contain the celebrated "ZePdovich approximation" 
(Zel'dovich 1970, 1973) as a special case, see Buchert 1989, 
1992) breaks down at the epoch of formation of caustics in 
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the density field. Fluid elements are treated as to propa- 
gate freely through each other; multi-stream flow arises 
beyond the caustic, i.e., the velocity field is no longer 
single-valued. Self-gravitation of the multi-stream sys- 
tem is the source for holding structures together: a con- 
siderable fraction of the particles is trapped within three- 
stream systems (called 'pancakes' in the cosmological con- 
text) which, in the course of time, develop a hierarchy of 
nested structures due to the formation of N-stream sys- 
tems with increasing N. Doroshkevich et al. (1980) have 
demonstrated this in a two-dimensional numerical simu- 
lation; the number of streams as a function of time was 
calculated by Kofman et al. (1994). This hierarchy of 
structures is also predicted by higher-order Lagrangian 
perturbation schemes (Buchert et al. 1997); compare also 
the works by Shukurov (1981) and Fillmore & Goldre- 
ich (1984a,b). The notion of 'non-dissipative gravitational 
turbulence' has been advanced for this phenomenon by 
Gurevich & Zybin (1995 and ref. therein). For further dis- 
cussions in the framework of Lagrangian perturbation the- 
ory on details of multiple shell-crossings related to bifur- 
cations of the Lagrangian solutions and the mathemati- 
cal problems involved see, e.g., Buchert & Ehlers (1993), 
Buchert (1994, 1996), Ehlers & Buchert (1997); in the rcl- 
ativistic context: Clarke & O'Donnell (1992). 

In order to repair this shortcoming, the "adhesion ap- 
proximation" (Gurbatov et al. 1989) has been proposed. 
This model proved to be successful, provided the struc- 
ture formation process is not followed too deeply into 
the nonlinear regime (Kofman et al. 1990, Weinberg & 
Gunn 1990a,b; Sahni & Coles 1995, Shandarin 1997). At 
late nonlinear stages the model predicts a distribution of 
large-scale structures which is different from that pre- 
dicted by N-body simulations, especially for much small- 
scale power. One reason for this was referred to the model 
itself, i.e. that the gravitational potential may acquire an 
additional part via nonlinear evolution on smaller scales 
(Kofman et al. 1992), it may also be traced back to the 
model's drawback of not conserving the linear momentum 
(Kofman & Shandarin 1988). We advocate the explana- 
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tion that it may also be a consequence of treating the 
"viscosity coefficient" as spatially constant, as will become 
clear later. A recent work by Shandarin & Sathyaprakash 
(1996) proposes a semi-numerical model which is built 
in the spirit of 'adhesion' of particles, but which has the 
advantage of conserving momentum by construction. In 
this paper we will construct models that allow follow- 
ing the structure formation process beyond shell-crossing 
time analytically having both the 'adhesive' property of 
'sticking' fluid elements and the conservation of momen- 
tum. 

We shall emphasize that the notion of a 'self-gravitating 
gas' (which we are going to specify) is more adequate when 
we speak about the regime of multi-stream flow in colli- 
sionless systems than the generally adopted simplification 
to "dust matter" : in the presence of multi-stream flow, 
one should focus on the average bulk motion rather than 
on individual particles. We shall argue that the action 
of the multi-stream system on the average motion can 
be modeled by means of pressure-like forces. Taking this 
point of view we will be able to propose an approximation 
scheme that is valid in the regime of multi-stream flow. 
The assumption of isotropic velocity dispersion singles out 
a model which determines a transport coefficient - not un- 
like the "viscosity" in the standard "adhesion approxima- 
tion" - which, however, is of gravitational origin. The ad- 
vantage of this point of view is two-fold: firstly, we do not 
have to invent a "viscosity term" to mimic gravitational 
adhesion of fluid elements after shell-crossing, but we can 
derive it from the gravitational action of the multi-stream 
system itself and, secondly, we can learn about the restric- 
tions we have to impose on the general problem in order 
to obtain models of this type. The latter will strongly sup- 
port our opinion that these restrictions should be relaxed. 
We wish to note that, starting from another reasoning, a 
similar interpretation of the "viscosity" in Burgers' equa- 
tion has been advocated by Ron Kates (priv. comm. and 
unpublished notes). 

We also think of several implications of such a descrip- 
tion, although we are not going to develop them in this pa- 
per. It is known that "dust" is only recovered as a singular 
limit: a pure single-stream system maintains a J-shaped 
velocity distribution along trajectories and generates a sin- 
gularity at shell-crossing. An initially small deviation from 
a ^-distribution allows, however, to follow the distribution 
function smoothly beyond shell-crossing; caustics occur- 
ring in the density field for "dust matter" are smoothed 
out by velocity dispersion (Zel'dovich & Shandarin 1982). 
Thus, besides the obvious advantages implied by a model 
which can be followed through shell-crossing, it also en- 
tails the description of how velocity dispersion evolves, 
which is a key element in the diagnostics of cosmic struc- 
ture. The "Jeans' length", which has not to be as small as 
the one associated with the usual thermostatic pressure, 
represents the scale on which structures are diffused (for 
a discussion of the Jeans' length in collisionless systems 



see Seliverstov 1974, for a simulation of neutrino distribu- 
tions see Melott 1982,1983, Ma & Bertschinger 1994), for a 
Fermi gas a lower limit is set by the phase space constraint 
(Trcmainc k Gunn 1979, Ruffini & Song 1987, Kull et al. 
1997). Including velocity dispersion will allow us to access 
smaller spatial scales than with previous analytical mod- 
els, and gives meaning to the notion of (virial) equilibrium 
in a collisionless system, which is a well-studied case in the 
context of stellar systems (Binney & Tremaine 1987). The 
treatment of multi-stream flow on cosmological scales will 
thus connect models for large-scale structure with small- 
scale systems which, apart from recent efforts to simu- 
late galaxy formation embedded into large-scale structure, 
have been treated as isolated systems so far. Yet another 
perspective is opened concerning reconstruction models 
of the density field from observed peculiar-velocity data 
(see Dekel 1994 for a review). Present reconstruction tech- 
niques suffer from the fact that they employ nonlinear 
models in regimes when singularities in the flow devel- 
oped; the quest for singularity-free models is especially 
obvious in that context (for a discussion of this issue see 
Susperregi & Buchert 1997). 

We proceed as follows: in Sect. 2 we introduce the sys- 
tem of equations that we are going to study, in Sect. 
3 we give an equation governing the evolution of the 
mean peculiar-velocity in the presence of multi-stream 
flow, discuss some properties of the dynamical evolution 
prescribed by this equation, and relate it to commonly 
employed models of large-scale structure. In Sect. 4 we 
present our criticism and perspectives. 

2. Basic equations in the presence of velocity 
dispersion. 

In this section we shall first use physical (i.e., non-rotating 
Eulcrian) coordinates. The introduction of so-called 'co- 
moving coordinates' (i.e., scaled Eulerian coordinates in- 
dexing fundamental observers in a background cosmology) 
will be employed when comparing with commonly used 
approximation schemes in cosmology in the next section. 
Let us consider a gas of N particles of mass m which in- 
teract only gravitationally. We may not only think of stars 
(which is the common point of view - see the textbooks 
by, e.g., Binney & Tremaine 1987, Ogorodnikov 1965 and 
Saslaw 1985 -), but alternatively of dark matter particles 
(e.g., non-relativistic massive neutrinos or WIMP's, e.g., 
Bertschinger 1993), or corresponding larger mass units 
(like galaxy halos) or, generally, of any patch of colli- 
sionless matter (a fluid clement or a Lagrangian particle). 
The interpretation of our model depends on the context to 
which we apply it. In any case we have to specify the na- 
ture of the particles to see whether the assumptions apply 
- at least in a phenomenological sense - (e.g., for a gas of 
galaxy halos we neglect merger events, or dynamical fric- 
tion; this also violates the particle number conservation; 
for any patch of matter we implicitly neglect additional 
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terms that would arise by coarse-graining the distribu- 
tion function, etc.). 

We define the one-particle distribution function /(x, v, t) 
such that /(x, v, t)dxdv is the probability density that 
a particle be within the volume element dx centered at 
the physical position x with a physical velocity within the 
volume element dv centered at v at time t. In the mean 
field approximation, it obeys Vlasov's equation (e.g., Bin- 
ney & Tremaine 1987) (summation over repeated indices 
is understood): 



dt % dxi ' % dvi 



(1) 



where the mean-field gravitational field strength g is de- 
termined self-consistently by the equations (the comma 
denotes spatial derivative, i.e., ^ = □§-): 

5M (x, t) = A - AnGmN f dvf(x, v, t) , (2a) 

9tj(x,t) =g jti (x,t) ; (2b) 

G denotes the gravitational constant and A the cosmolog- 
ical constant. 

We now define the averaging operation (...}: If 0(v) is 
any function of the velocity v, then 



(0):= 



mN 
g(x,t) 



dv<l)(v)f(x,v,t) 



(3) 



where the mass density is g — mN fdvf. We define 
the average velocity v := (v), and the stress tensor 
II = £>((vv) — vv), which takes velocity dispersion into 
account, i.e., the fact that at a given position there are 
several particles with different velocities. (It should be no- 
ticed that II is a positive definite, symmetric tensor). 

From Vlasov's equation (|l|) the equations obeyed by 
the density g, the average velocity v and the stress tensor 
II are found by integrating out (||), Uj-(0), and ^^-(0) 
over velocity space: 



velocity and so on. We thus get an infinite hierarchy of 
equations for all moments of the velocity. 

The task is now to find a way to close this hierarchy. 
In hydrodynamics the common practice is to assume that 
collisions lead to local equilibrium on a short time-scale, 
so that the hydrodynamical system of equations is reduced 
to balance equations on a much larger time-scale for the 
five collisional invariants g, Vi and the internal energy U . 
In dilute gases, the Chapmann-Enskog expansion is then 
applied to find Navier-Stokes' equation and Fourier's law 
(e.g., Balescu 1991). This method cannot be employed in 
our case, because collisions are neglected and local equi- 
librium is thus undefined. 

The simplest way to close the hierarchy without re- 
sorting to local equilibrium considerations is to neglect 
velocity dispersion altogether, i.e., / oc 6(v — v) which is 
the definition of "dust matter" (e.g., Gurevich & Zybin 
1995). Then By = 0, L ijk = and Eqs. (4) reduce to 



d t g+ (QVi) yi = 
d t Vi +v j v ii j = g t 



(5a) 
(5b) 



which is the well-known set of equations (together with 
(2)) employed to model gravitational structure formation, 
being equivalent to following the motion of individual par- 
ticles. 

As remarked in the introduction, these equations have 
the problem that their solutions develop caustics, precisely 
because of neglecting velocity dispersion. We go one step 
further by explicitly including velocity dispersion from the 
beginning. Let us therefore assume that there is a small 
velocity dispersion of order e << 1. This means that the 
typical relative deviations of particle velocities from the 
average velocity v are of order e (i.e., \vt yp — v| ~ e|v|). 
This assumption implies that Hij ~ ££ 2 |v| 2 and Lijk ~ 
g£ 3 |v| 3 . To study the effect of velocity dispersion, we shall 
retain only the leading order terms due to it, i.e. we shall 
work only up to order e 2 . [J Then Lijk : k can be dropped 
in Eq. ( |4c[ ) provided Lijk changes spatially on the same 
scale as Vi and we get (together with (2)) a closed set of 
equations for g, v, B and g: 



d t g+ (QVi) t i = , 
d t Vi + VjVij =gi Tlijj , 



(4a) 
(4b) 



d t Hij + VkHij.k + Vk,k^hj — ~TlikVj,k — KjkVi,k — Lijk.k , 

(4c) 



where L ijk = g((vi - «»)(«,- - Vj)(vk - «fc)}- Eq. Q) is 
known as Jeans' equation in stellar systems theory. This is 
not a closed set of equations for the density g, the average 
velocity v and the stress tensor B and (together with (2)) 
the gravitational field strength g, because oiLijk- We may 
in fact find the equation obeyed by this third-rank tensor, 
but it turns out to depend on the fourth moment of the 



d t g+ (gv l ) : , = 
d t Vi +VjVi > j = g 



in, 



9*B 



■f-Hj 



Vk^-ij,k + Vk,k^-ij — —H-ikVj : k — BjfeW^fc 



(6a) 
(6b) 
(6c) 



Using this set of equations we go one step further down 
the hierarchy of equations for the velocity moments. It is, 
however, still difficult to cope with these equations and, in 
the present paper, we introduce the strong simplification 



It must be remarked that this is not equivalent to assum- 
ing that the velocities be Gaussian distributed. Gaussianity 
imposes an additional constraint on the set of Eqs. (6). 
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that the velocity dispersion is approximately isotropic, i.e. 
the stress tensor 11^ is diagonal and has only one indepen- 
dent component (a pressure-like term): 



Combining this equation with the continuity equation (6a ) 
we find 



n, 



' ph. 



p > 



(7) 



This assumption is certainly sensible at early stages 
of the structure formation process: we can think of an 
initial condition with isotropic velocity dispersion (it has 
to be so in the case of a homogeneous-isotropic matter 
distribution). Then, at earlier stages of the evolution, it 
will remain approximately isotropic. At later stages we 
have to consider the isotropy of the stress tensor as an 
idealization of the generally anisotropic velocity ellipsoid. 
This assumption is to be considered a weak point in our 
analysis, since there is nothing to prevent the velocity dis- 
tribution from becoming largely anisotropic in contrast 
to the situation in a collisional system; it may work for 
single galaxy halos, but the bulk of the medium is not 
isotropic^]. Isotropy may, however, provide a good model 
if violent relaxation is taken into account (see Lynden-Bell 
1967, Henriksen & Widrow 1997, White 1996 and Kull et 
al. 1997). 
With this assumption, Eq. ( pq ) is greatly simplified: 



P + -xVi yi p 



2paij = 



(8) 



where we have introduced the total or Lagrangian deriva- 
tive along integral curves of the mean velocity field, 
:= dt + Vidi, and the shear tensor o~ij by splitting the 
symmetric part of the mean-velocity gradient: Vuj) — 
\{Vij + Vj,i) = %Vk,k$ij + o- t j. 

Notice that Eq. m) imposes a strong kinematical re- 
striction, namely, o~ij = 0. Although we generalize the 
basic system of equations by allowing for a nonvanishing 
p, the interpretation of p in terms of a strictly isotropic 
velocity dispersion will accordingly restrict the kinemat- 
ics of the mean motion. This makes clear that, unless we 
do not consider p as a phenomenological generalization of 
the "dust matter model", the assumption (J7|) is far too 
restrictive. It is, however, one of the assumptions which 
has to be imposed in order to recover the "adhesion ap- 
proximation" . 

Keeping this in mind we go on by evaluating the trace 
of Eq. (||); it yields a relation between the density and the 
"dynamical pressure" : 



p = d t p + Vip. 



-,PVi 



(9) 



2 From the observational point of view velocity dispersion 
can only be inferred from the line-of-sight information, the 
isotropy assumption is then employed when concluding on the 
3D dispersion value. 



5p . 

P= -X-Q 
3q 



(10) 



For this equation it is straightforward to find its gen- 
eral (Lagrangian) integral (which we may interprete as 
an 'equation of state' for a single volume element): 



p(X,t) = K(X) e (X,t) 



5/3 



(11) 



where k is positive and of order e 2 , and X are Lagrangian 
coordinates which label fluid elements and coincide with 
the Eulerian coordinates x at the initial time t = to. (Note 
that the integration of ( |1C| ) implies that k^O; k = is a 
singular limit of the equations under consideration.) 

If the physical system were a true fluid, (|IJ) would be 
interpreted as describing the adiabatic evolution of each 
volume element in an ideal gas (consistently with the fact 
that we have assumed T-njh,k — 0, which would be in- 
terpreted as the absence of heat flux between neighbour- 
ing volume elements). But unlike in a fluid, there are no 
collisions in the system we are studying that prevent ve- 
locity dispersion from growing and/or becoming highly 
anisotropic. 

In the sequel we restrict the general integral (O) to the 
case where the initial data p(X, to) and p(X, t ) are chosen 
such that k is independent of X. Therefore, we drop the 
freedom of givingp(X, t ) independently of g(X, t ). Thus, 
along any trajectory, we have the relationship p — Kg 5 ' 3, 
with the same constant k. Consequently, we have the same 
relation also at any point in Eulerian space. The integral 
then attains the status of a (global) 'equation of state' for 
the medium occupying Eulerian space. However, here, it 
is not imposed as an equation of state in thermodynamics, 
but follows from the dynamical equations under consider- 
ation and the assumptions employed. 



Taking this result into Eq. (|6b|) we get our final set of 
equations: 





d t g + (gvi) 


4=0 , 


(12a) 




d t Vi +VjVi,j = gi - 


5 —l /a 
- x^g ' g,i , 


(12b) 


9i,i 


= A - AiiGg , 


9i,j = 9j,i ■ 


(12c) 



The difference between the closed set of Eqs. (12) and the 
one usually employed in cosmology (2,5) is the pressure- 
like term oc Q~ l ' 3 Q,i, which takes into account (under the 
approximations stated above) velocity dispersion (i.e., the 
dynamical reaction due to multi-stream regions). 
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3. Equation for the mean peculiar velocity 
relative to a background Hubble flow. 

We shall now show what novel features arise in the dy- 
namical evolution prescribed by (12) due to the inclusion 
of the pressure-like term in contrast to the evolution pre- 
scribed by (2,5). 

In order to compare our model with the commonly 
studied approximation schemes for large-scale structure 
we shall perform a change of variables. We consider the 
homogeneous-isotropic solutions of the basic equations 
(Friedmann-Lemaitre backgrounds) characterized by the 
expansion factor a(t) (and Hubble's function H = a/ a), 
and define 'comoving coordinates' q, an average peculiar- 
velocity field u and a (mean field) gravitational peculiar- 
acceleration w as follows: 

q:=-x,u:=v-iJx,w:=gH x, (13) 

a 3 

where Qn(t) is the homogeneous background density (sat- 
isfying Qh + 'SHqh = 0), which coincides with the spatially 
averaged density, if we impose periodic boundary condi- 
tions on u and w; in this case the Hubble-flow exists and 
is uniquely defined (see Buchert & Ehlers 1997)Q 

In terms of these variables, Eqs. (12) become (spatial 
derivatives refer now to q and time derivatives are taken 
at constant q; hereafter, we drop the bar above u for no- 
tational simplicity): 



d t Q + 3HQ+ -{eui),i = , 
a 

r. 1 TT 5K _,/, 

d t Ui + -UjUij +Mui = Wi - --q ' g.i , 



Wi, 



-AirGa(g - g H ) 



Wi. 



Wi 



(14a) 

(14b) 
(14c) 



Combining ( |l4q ) and (14c) we find the following equa- 
tion, using the Lagrangian derivative operator ' := d t \ x + 
Vidi\ x = dt\ q + — di\q (written in vector form): 



u + Hvl = w + (Aw , 

-1/3 



(15) 



with the coefficient £ = § 4^g^r g~ ' > 0, which depends 
on density and explicitly on time. The difference between 
Eq. (|lj) and the one generally used to model large-scale 
structure formation (which is found by setting £ = 0) is 
the Aw term. 

Recall first the case of "dust" (i.e., no velocity disper- 
sion, C = 0). In the weakly nonlinear regime Zel'dovich's 
approximation (Zel'dovich 1970, 1973) is a successful 

3 It is straightforward to show that this result carries over to 
the case where velocity dispersion is present, since the Hubble 
flow cancels in the expression for II, II = g((w) — vv) = 
p((uu) — uu), i.e., II is a tensor field on the torus, and the 
divergence of the tensor appears in the basic equations, i.e. the 
corresponding flux vanishes by integrating over the torus. 



model until shell-crossing singularities develop. The tra- 
jectories in that approximation obey the parallelism of 
peculiar-gravitational acceleration and peculiar-velocity 
(see, e.g. Bildhauer & Buchert 1991, Kofman 1991, 
Buchert 1992), 



w = F(t)u , F{t) = 4nGg H ^ 



b(t) 



b(t) 



(16) 



where b(t) is the growing mode solution of the linear the- 
ory of gravitational instability for "dust" (i.e., it solves 
the equation b + 2Hb — AnGgnb = 0). 
Since Eqs. (14) were obtained under the condition of small 
velocity dispersion, we can try to extrapolate Zel'dovich's 
approximation (|l6[) into this regime (which is equivalent 
to solving Eq. ( [15|) by iteration) and thus find from Eq. 



u+(H -F)u = (F{t)Au 



(17) 



In order to construct the model we have to derive from 
the solution of Eq. ( jl7| ) the trajectory field of the average 
flow q = F(X, t) by quadrature: u = a¥. Changing the 
temporal variable from t to b (this is possible since b(t) is 
a monotonically increasing function of time) and defining 
a rescaled velocity field u = u/ab, Eq. (|l7|) becomes 



eiu 



^ db ''~ db 



(18) 



where \i = QF{t)/b > 0. If /i were independent of den- 
sity, this equation would become 3D Burgers equation, 
whose solution is analytically known. The fact that /i de- 
pends on density presents an obstacle for finding an ana- 
lytical solution. The principal advantage of Eq. ( |l8| ) over 
the same equation with (i = is that it does not lead to 
caustic formation, since velocity dispersion smoothes out 
the singularity (Zel'dovich & Shandarin 1982, Shandarin 
& Zel'dovich 1989; see Ginanneschi 1998 and ref therein 
for a thorough analysis). Therefore, this equation could 
be used as it stands to follow the dynamical evolution be- 
yond the time when singularities in the "dust" continuum 
would arise. To be more precise, "shell-crossing" would 
still happen, but this doesn't manifest itself as a singular- 
ity in the average flow, rather as a (smooth) peak in the 
density field. Morphologically distinct patterns, classified 
by the Lagrange-singularity theory in the case of "dust" 
(see Arnol'd et al. 1982), will also emerge in the sense of 
smoothed-out images of critical sets on the Lagrangian 
manifold of the "dust" medium. This picture might not 
be true for large velocity dispersion. 

As time goes by and the system becomes more and 
more virialized, Eqs. (14) cease to be a good description 
of the dynamical evolution, because velocity dispersion 
generically both grows and becomes anisotropic. 
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Eq. (|18j) is formally equivalent to the key equation of 
the "adhesion model" in which the coefficient \i is positive 
and constant in space and time (Gurbatov et al. 1989), and 
we see that the Au term behaves as a viscous force. In fact, 
in the singular limit /i — > (the so-called inviscid limit), 
we approach a singular case of the "adhesion model" (Gur- 
batov et al. 1989; see also Gurbatov et al. 1983, 1985), 
the so-called "sticky particle model" for which geometri- 
cal construction techniques have been advanced (Pogosyan 
1989, Kofman et al. 1990, Sahni & Coles 1995). For further 
details on models based on Burgers' equation the reader 
may also consult the book by Gurbatov et al. (1991) and 
the interesting article by Vergassola et al. (1994). Further 
insight into possible applications is provided by a model of 
Jones (1996) for a two-component collisionless-baryonic 
system, which is related to the "adhesion model" . 

Although we here recover a variant of the "adhesion 
model", it should be clear that the new model presented 
implies an improvement for various reasons. The main one 
is that the dissipative-like term appears in a natural way 
together with a clear physical interpretation (i.e., velocity 
dispersion) , unlike in the "adhesion model" , where it is 
motivated by phenomenological arguments. 

This fact also solves the 'momentum-conservation vi- 
olation problem' arising in the "adhesion model" (Kof- 
man & Shandarin 1988; see also Gurbatov et al. 1983, 
Shandarin & Sathyaprakash 1996). Indeed, if velocity dis- 
persion is neglected, Eq. ( J5b[ ) describes the evolution of 
the velocity of individual particles, so that appending a 
Au force to this equation implies that momentum is not 
conserved. On the other hand, the velocity field u(q, t) is 
the mean velocity of the particles at the (comoving) Eu- 
lerian position q at time t. For the mean velocity, gravity 
is not the only force, there is also an effective force due 
to velocity dispersion, so that the momentum balance is 
not violated. This effective force (pressure) simply pro- 
vokes a flow of kinetic energy between the bulk average 
motion and the "random" motion of particles (described 
by velocity dispersion). The total kinetic energy (the ki- 
netic energy of the mean flow plus internal kinetic energy) 
changes only due to gravitational work. 

The conservation of mass and momentum is already 



evident from the way we have found Eqs. (14a) and (jT4 



which are balance equations expressing the conservation 
of mass and momentum. In fact, since the pressure force 
appears as a gradient (in general, as the divergence of 
a tensor for anisotropic velocity dispersion), it describes 
the transfer of momentum between neighbouring volume 
elements. 

The reversibility of the transfer of momentum can also 
be immediately appreciated from our model equations. 
The coefficient /i is not related to an irreversible process: 
if one performs a time reversion t — > —t, u — > — u (w, g 
remain unchanged), then F — ► — F, £ — > £, and Eq. (W% 
is invariant under time reversion. (This is true also for the 
exact equation (Ha) as well as for the basic system of equa- 



tions. Therefore, the coefficient £ should not be considered 
as a viscosity, a true viscosity coefficient would not change 
sign under time-reversion. Another difference is that the 
coefficient \x depends on the initial conditions (through n) , 
whereas true damping rates are insensitive to initial con- 
ditions. We may call the coefficient /i gravitational multi- 
stream or GM-coefficient, because it arises from the self- 
gravitation of multi-stream systems. The phenomenon of 
"reversible damping" is not unknown in physics (e.g. Lan- 
dau damping); the phenomenon we are discussing coin- 
cides with 'non-dissipative gravitational turbulence' as de- 
scribed by Gurevich & Zybin (1995). 

Another interesting point is that the GM-coefficient /x 



in (18) is not arbitrary: it has a dependence on time and 



density which is determined by the dynamical equations. 
However, this (very plausible) fact complicates the model 
which no longer can be solved by the known solution of 
the 3D Burgers equation. We have not found an immedi- 
ate variable change to perform Cole-Hopf-type or other 
transformations in order to reduce the model to the lin- 
ear diffusion equation. Application of the methods which 
work in the case of the 3D Burgers equation to more gen- 
eral problems apparently creates obstacles (see, e.g., Ner- 
ney et al. 1996). The model must be solved numerically 
by determining /i(g>, t) at each time-step for each particle. 

4. Discussion and criticism. 

We have presented a system of equations that allows 
studying the structure formation process at stages where 
the usually employed Euler-Poisson system breaks down. 
We discussed the possibility of constructing models for 
large-scale structure which do not suffer from the occur- 
rence of shell-crossing singularities. We also proposed a 
model which features "adhesion" of fluid elements sim- 
ilar to the commonly used "adhesion approximation", 
and allows following the structure formation process be- 
yond stages where a laminar fluid approximation based 
on "dust" breaks down due to the development of multi- 
streaming: a "viscous" term appears in the model equa- 
tion which we derived from the combined action of multi- 
stream flow and self-gravity. 

Besides the obvious advantages that we do not have to 
invent an adhesive term and that momentum conservation 
is not violated, we also see a further merit of this approach: 
we have gained some physical insight into the "adhesion 
model" and could even improve on it. 

One may ask why the combined effect of self-gravity 
and velocity dispersion actually "looks like" a dissipative 
term. To clarify this, let us add the following illustration: 
consider a vessel, filled with a non-viscous toy fluid, and 
placed into the gravitational field of the Earth. Initially, 
the fluid may be at rest and its density homogeneous. 
This state is clearly unstable and the fluid will evolve into 
a stable, inhomogeneous and stratified state. The overall 
kinetic energy vanishes in this example, but the gravita- 
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tional potential energy has decreased. The "missing en- 
ergy" has not been dissipated away, rather it has been 
transformed into internal kinetic energy due to compres- 
sion work; the final volume of the fluid is smaller than the 
initial volume, and the pressure has always opposed the 
motion induced by gravity. This is similar to what happens 
in the cosmological setting: volume elements would have 
equal tendency to compress or expand due to pressure- 
like forces, but the presence of gravity favors contraction 
of volume elements, especially when they are falling onto 
high-density regions. The effect of transformation of en- 
ergy into internal velocity dispersion is large, because not 
only potential energy, but also bulk kinetic energy de- 
creases, when the fluid elements crash onto high-density 
walls. 

We have also shown that pressure-like forces may be 
relevant even on cosmological scales and an effective Jeans' 
length may be larger than the (comoving) globular cluster 
scale, which is commonly attributed to a hydrodynami- 
cal (ideal gas) pressure (compare this to the effect of a 
gas pressure within pancakes investigated by Sunyaev & 
Zel'dovich 1972). However, we now understand clearly the 
limited status of this and similar models. We shall now 
elaborate on this criticism in more detail. 

The extrapolation of the relationship (|l^) into the 
weakly non-linear regime is well justified for "dust". If 
we define the weakly non-linear regime as suggested in 
(Buchert 1989), namely, as a linearization in Lagrangian 
space, Zel'dovich's extrapolation of the Eulerian linear 
perturbation solution into the nonlinear regime is not only 
consistent with Lagrangian linearization of the full Euler- 
Poisson system (Buchert 1989, see also Doroshkevich et 
al. 1973 for a self-consistency test for part of the system), 
but arises naturally as a subcase of first-order solutions 
in a Lagrangian perturbation approach (Buchert 1992). 
This self-consistency of the "dust model" is mirrored in 
the fact that the trajectories of Lagrangian first-order per- 
turbations even provide a class of 3D exact solutions to 
the Euler-Poisson system (Buchert 1989). 

If pressure is taken into account, extrapolation of the 
"dust" model trajectories (as we had to do in order to 
recover the "adhesion approximation" ) is no longer per- 
mitted. As we have shown it may be used as a first ap- 
proximation for small velocity dispersion (and this is pre- 
cisely the method followed in Zel'dovich & Shandarin 1982 
and Bharadwaj 1996). A similar remark applies to other 
assumptions of extrapolation like the Frozen-Potential- 
Approximation (Brainerd et al. 1993, Bagla & Padman- 
abhan 1994), which we may employ instead of (plf). This 
approach simply follows the evolution of velocity disper- 
sion along the trajectories dictated by the "dust" model, 
its main drawback being that the "back-reaction" of ve- 
locity dispersion on these trajectories is not taken into 
account. In fact, the simple version of parallelism (11(3) is 
not even justified in the (Eulerian) linear regime (Buchert 
et al. 1998). The break-down of the assumption ( |16[) after 



pancake formation is plausible and it is well-known for 
epochs after shell-crossing to which we want to apply the 
model (Doroshkevich 1973). 

Nevertheless, the "adhesion approximation" works rather 
well, because the model is taken at face value and not 
restricted a posteriori by the constraints which led us to 
its derivation. Accordingly, we expect the model we pro- 
posed to improve on the "adhesion approximation" , since 
it specifies the dependence of "adhesion" on the local den- 
sity. The publication of this work was motivated by this 
gain of insight and improvement. Let us summarize the 
constraining assumptions which restrict the general prob- 
lem to the "adhesion approximation" : 

• small velocity dispersion; we keep only the first equation 
in the hierarchy of velocity moments (Al). 

• isotropy of the dispersion tensor which implies isotropy 
of the mean motion (A2). 

• parallelity of peculiar-velocity and -acceleration, both 
defined relative to a global Hubble-flow (A3). 

• a spatially constant relationship between initial density 
and initial pressure (A4). 

• the GM-coefficient is constant in space and time as well 
as positive (A5) (which could be achieved by imposing a 
relationship p = Kg 2 ). 

This list of assumptions shows that there is enough room 
for a generalization of the "adhesion approximation" . 
Our main concern should be focussed on assumptions (Al) 
and (A2). Although both of them may be considered a 
good working hypothesis from a theoretical point of view, 
we know that the physical situation is not in favor of these 
assumptions, especially if we follow the average flow fur- 
ther into the nonlinear regime. Velocity dispersion will 
become quickly large in the multi-stream regime. A hint 
of this can be seen in observations of velocity dispersion in 
rich clusters of galaxies which is of the order of 1000 km/s, 
i.e. certainly not small compared with the bulk speed. 
However, viewing assumptions (Al) and (A2) together, we 
would obtain substantially similar models by arguing on 
phenomenological grounds: assumption (A2) gives an ide- 
alized model for the dispersion tensor and it just needs a 
relationship p = j3{g) to close the hierarchy of velocity mo- 
ments. In light of a phenomenological approach assump- 
tion (Al) or, alternatively, assumption (A5) just specify 
the function (3, but we may as well consider other relation- 
ships as discussed in (Buchert et al. 1998). This problem 
calls for a better model that relates the dispersion tensor 
to the density and velocity fields, but it does not invalidate 
the main ideas exemplified with assumption (Al). 

Concerning assumption (A3) a first indication of how 
to proceed in order to investigate a valid extrapolation 



into the weakly nonlinear regime is furnished by Eq. (15) 



itself: the adhesive term is proportional to Aw rather than 
Au. It is possible to derive an evolution equation for the 
peculiar-gravitational field strength w. Its specification to 
the weakly nonlinear regime, however, requires Lagrangian 
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perturbation techniques that lie beyond the scope of the 
present work. Assumptions (A4), (A5) have their origin in 
mathematical simplicity and could be easily relaxed in a 
numerical simulation of the model. 
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